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Abstract 

Background: The mesic habitats of eastern Australia harbour a highly diverse fauna. We examined the impact of 
climatic oscillations and recognised biogeographic barriers on the evolutionary history of the delicate skink 
{Lampropholis delicata), a species that occurs in moist habitats throughout eastern Australia. The delicate skink is a 
common and widespread species whose distribution spans 26° of latitude and nine major biogeographic barriers 
in eastern Australia. Sequence data were obtained from four mitochondrial genes (ND2, ND4, 12SrRNA, 16SrRNA) for 
238 individuals from 120 populations across the entire native distribution of the species. The evolutionary history 
and diversification of the delicate skink was investigated using a range of phylogenetic (Maximum Likelihood, 
Bayesian) and phylogeographic analyses (genetic diversity, <I> ST , AMOVA, Tajima's D, Fu's F statistic). 

Results: Nine geographically structured, genetically divergent clades were identified within the delicate skink. The 
main clades diverged during the late Miocene-Pliocene, coinciding with the decline and fragmentation of 
rainforest and other wet forest habitats in eastern Australia. Most of the phylogeographic breaks within the delicate 
skink were concordant with dry habitat or high elevation barriers, including several recognised biogeographic 
barriers in eastern Australia (Burdekin Gap, St Lawrence Gap, McPherson Range, Hunter Valley, southern New South 
Wales). Genetically divergent populations were also located in high elevation topographic isolates inland from the 
main range of L delicata (Kroombit Tops, Blackdown Tablelands, Coolah Tops). The species colonised South 
Australia from southern New South Wales via an inland route, possibly along the Murray River system. There is 
evidence for recent expansion of the species range across eastern Victoria and into Tasmania, via the Bassian 
Isthmus, during the late Pleistocene. 

Conclusions: The delicate skink is a single widespread, but genetically variable, species. This study provides the 
first detailed phylogeographic investigation of a widespread species whose distribution spans virtually all of the 
major biogeographic barriers in eastern Australia. 



Background 

The coastal regions of eastern Australia are currently 
dominated by wet forest and drier sclerophyllous habi- 
tats that harbour a highly diverse fauna [1,2]. While the 
majority (-70%) of the Australian continent is covered 
by arid or semi-arid vegetation, eastern Australia pro- 
vides a narrow, but largely continuous expanse of habi- 
tat for mesic-adapted species [1,3,4], These mesic 
habitats are generated through the presence of the 
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Great Dividing Range (GDR), which abuts the entire 
length of the east coast (-2,500 km) in a north-south 
alignment ([5-7]; Figure 1). In the context of an expan- 
sive continent that is characterised by low topographic 
relief, the moderate elevation (-1000-1300 m, maximum 
-2300 m) provided by the GDR generates altitudinal, 
climatic and environmental variation, and precipitates 
the required moisture to support mesic vegetation 
[3,5,7], 

Although widespread glaciation never occurred in 
Australia [8,9], climatic oscillations have driven repeated 
altitudinal and distributional shifts in mesic habitats 
along the eastern margin of the continent (reviewed in 
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Figure 1 The major biogeographic barriers in eastern Australia. A description of each barrier is provided in Table 1. The location of 
Kroombit Tops is shown in Figure 2. Inset: The coastline of the Bass Strait region 14 kya. Tasmania has repeatedly been connected to the 
mainland during the Pleistocene by land bridges, with the most recent connection occurring 43-14 kya during the last glacial maxima. The 
western land bridge was severed 17.5 kya, with the eastern connection (the Bassian Isthmus) being inundated 13 kya, isolating Tasmania from 
the mainland (after [16]). 
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[4]). Palaeoclimatic studies indicate that the extent and 
composition of the vegetation has fluctuated dramati- 
cally over the last 10 myr, although there has been a 
general transition from rainforest towards drier environ- 
ments and sclerophyllous vegetation [1,3,8]. The rainfor- 
ests that had previously dominated eastern Australia 
contracted between the mid- and late-Miocene, giving 
way to woodland and open forest vegetation that was 
more suited to the drier climates [3,10-12]. Lowered sea 
level associated with globally drier conditions facilitated 
the expansion of vegetation into the low lying regions of 
south-eastern Australia (e.g. Gippsland Basin, Murray 
Basin; Figure 1) that had previously been subject to mar- 
ine inundation [3,11,12]. Although the extent of rainfor- 
ests briefly expanded again during the early Pliocene 
due to a temporary return to warm and wet conditions, 
by the end of the Pliocene open woodlands, sclerophyl- 
lous forests and grasslands dominated the landscape of 
eastern Australia [3,6,8,10]. 

The cool-dry to warm-wet climatic fluctuations that 
commenced during the Pliocene intensified throughout 
the Pleistocene and led to the repeated expansion and 
contraction of mesic habitats in eastern Australia and 
the regular encroachment of drier habitats into the 
coastal fringes [4,8,13,14]. There was periodic flooding 
of the low lying coastal and inland basins in eastern 
Australia during the sea level changes associated with 
these climatic cycles [3,6,15], which also resulted in the 
connection of Tasmania (TAS) to the mainland during 
glacial periods by Bass Strait land bridges [16] (Figure 
1). At present, the once widespread rainforest and wet 



forest vegetation is restricted to small, scattered rem- 
nants within a mosaic of dry sclerophyll woodlands and 
open forests along the east coast [1,17]. 

The evolutionary history of the resident fauna of the 
narrow mesic strip along the east coast has been influ- 
enced by both habitat barriers and physical barriers (e.g. 
mountain ranges, sea straits), which led to genetic diver- 
gence and, in some cases, speciation of allopatric popu- 
lations [1,18]. The most well-studied barrier in eastern 
Australia has been the Black Mountain Corridor (BMC) 
in the Wet Tropics region of north Queensland (QLD). 
This thin strip of rainforest currently connects the 
northern and southern rainforest block of the Wet Tro- 
pics but was repeatedly severed in the past by dry forest 
habitats during globally drier climates [18,19]. Intensive 
research has revealed largely concordant patterns of 
genetic divergence across the barrier in a wide range of 
rainforest taxa [e.g. [19-24]], and improved our under- 
standing of how these barriers, in concert with climatic 
oscillations, have generated the high levels of biodiver- 
sity evident in eastern Australia [18,25]. However, at 
least nine other biogeographic barriers have been identi- 
fied in eastern Australia (Tables 1 and 2; Figure 1), sev- 
eral of which have yet to be investigated in detail. These 
include dry habitat barriers (Burdekin Gap, St Lawrence 
Gap, Hunter Valley), mountain ranges that act as topo- 
graphic barriers (McPherson Range, southern New 
South Wales [NSW]), disjunct inland mountains 
(Kroombit Tops), sea straits (Bass Strait), and marine 
basins (Gippsland Basin, Murray Basin) (Tables 1 and 2; 
Figure 1). 



Table 1 Description of the recognised biogeographic barriers in eastern Australia (see Figure 1) 



Barrier 



Explanation of Barrier 



Burdekin Gap A broad region of dry woodland and savanna that extends to the coast and delineates the boundary between the northern 
rainforests and the mid-eastern Queensland forests [93-95] 

St Lawrence A dry habitat corridor that separates the mid-eastern Queensland forests from the south-eastern Queensland forests [93-95] 
Gap 

Kroombit Tops A disjunct inland region of high elevation moist habitat that is surrounded by drier eucalypt woodland. An inland cool and wet 
refuge for rainforest and wet forest adapted species [70,71] 

McPherson An east-west spur of the predominately north-south Great Dividing Range that runs along the Queensland/New South Wales 
Range border. A montane block of wet forest that represents a hybrid zone for birds and a barrier for lowland and dry forest plant 

species [93-96] 

Hunter Valley A dry, open, lowland river valley that delineates the southern limit of the eastern biogeographic region and the northern limit of 
the south-east forest region [93-96] 

Southern NSW Transition from the lowland coastal region to the higher elevation southern highlands region of the Great Dividing Range in New 
South Wales [71,74] 

East Gippsland Low lying coastal region that has been subject to repeated marine incursion (i.e. Gippsland Basin); abutted to the north by higher 
elevation regions of the Great Dividing Range [16,97] 

Bass Strait The shallow sea strait (depth 50-80 m, width 240 km) that separates Tasmania from mainland Australia. Land bridges have 

periodically connected the two landmasses during Pleistocene glacial periods (see Figure 1), with the last connection severed 13 
kya [16] 

Murray Basin Low lying region that has been subject to repeated marine incursion (i.e. Murray Basin), bordered to the west by the Mt Lofty 
Ranges, a known refugia [92,93,97] 

The Black Mountain Corridor in the Wet Tropics of north Queensland is not included as the distribution of the delicate skink (Lampropholis delicata) does not 
span this biogeographic barrier. 
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Table 2 The impact of the recognised biogeographic barriers in eastern Australia (Figure 1, Table 1) on vertebrates, 
invertebrates and plants 

Taxa Burdekin St Lawrence Kroombit McPherson Hunter Southern East Bass Murray 

Gap Gap Tops Range Valley NSW Gippsland Strait Basin 

Vertebrates 

Freshwater fish 

Pseudomugil signifier YN - S N - - - - 

[69] 

Amphibians 



Crinia signifera [79] 






s 


S 


N 


Y 


Y 


Y 


Y 


Limnodynastes peronii 
[67] 


Y 


N 


s 


N 


Y 


S 


S 


N 




Limnodynastes 
tasmaniensis [67] 


Y 


N 


s 


N 


S 


S 


Y 


N 


Y 


Litoria aurea [98] 


- 


- 


- 


- 


N 


- 


- 


- 


- 


Litoria citropa species 
group [73,77] 


- 


- 


Y 


Y 


Y 


Y 


- 


- 


- 


Litoria fallax [66] 


Y 


S 


Y 


Y 


- 


- 


- 


- 


- 


eptiles 




















Acritoscincus duperreyi 
[76] 














S 


Y 


Y 


Acritoscincus 
platynotum [76] 


- 


- 


- 


- 


Y 


S 


- 


- 


- 


Carlia rubrigularis/ 
rhomboidalis [22] 


Y 


- 


- 


- 


- 


- 


- 


- 


- 


Diporiphora australis 

|OuJ 


Y 


Y 


S 


N 


- 


- 


- 


- 


- 


Hoplocephalus 
stephensi [78] 






s 


Y 












Lampropholis 
gulchenotl [26] 








v 
Y 


v 
Y 


v 

T 


Y 




v 
Y 


Lerista bougainvilii [83] 










s 


S 


s 


Y 


N 


Liopholis whitii [74] 








s 


Y 


Y 


Y 


K 1 


Y 


Notechis scutatus [85] 








s 


S 


S 


S 


K 1 


Y 


Saproscincus 
mustelinus [75] 










Y 


S 


s 






irds 




















Ptilonorhynchus 
violaceus [71] 


S 


S 


Y 


N 


N 


Y 


N 






Sericornis frontalis [1 9] 


N 


N 


S 


S 


S 


S 


S 




S 


Sericornis magnirostris 
[19] 


Y 


M 


S 


S 


S 










Sericornis citeogularis 
[19] 


Y 




S 


S 


S 










lammals 




















Dasyurus maculates 
[84] 


N 


N 




s 


N 


s 


s 


Y 




Petaurus australis [99] 


Y 


S 


S 


N 


S 


s 


s 







Invertebrates 

Catomerus polymerus Y 
(M) [88] 

Catostylus mosaicus Y 
(M) [86] 

Drosophila birchii [1 00] N S 

Nerita atramentosa (M) Y 
[87] 

Plants 

Eucalyptus grandis S S - N S 

[101] 

Impact codes: Y = genetic break present across barrier; N: no genetic break observed across the barrier; S = insufficient sampling to examine the impact of the 
barrier;-= species distribution does not span the barrier. Several marine species (M) are included to investigate the impact of Bass Strait land bridges. 
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A recent study investigated the impact of five biogeo- 
graphic barriers in south-eastern Australia [26]; how- 
ever, here we adopt a broader approach and examine 
the influence of nine biogeographic barriers (Tables 1 
and 2) throughout eastern Australia on the evolutionary 
history of the resident biota. In particular, we focus on 
the delicate skink, Lampropholis delicata (De Vis, 1888 
[27]), which is unusual in that its distribution is so 
broad that it spans all of these barriers in eastern Aus- 
tralia (Figure 2, 3). The delicate skink is a small lizard 
(adult snout-vent length 35-51 mm) whose distribution 
extends across 26° of latitude from Cairns in north QLD 
to Hobart in TAS, with disjunct populations in far wes- 
tern Victoria (VIC) and south-eastern South Australia 
(SA) ([28]; Figure 2, 3). It is a common species that 
occurs across a range of moist habitats, including rain- 
forests, wet sclerophyll forests, woodland and heaths 
[28]. However, it also thrives in disturbed habitats and is 
one the most common skink species in suburban gar- 
dens along the east coast [28]. 

Here we examine the phylogeography of the delicate 
skink using 2426 bp of mitochondrial DNA sequence 
data (ND2, ND4, 12SrRNA, 16SrRNA) from across the 
entire native range of the species (Figure 2, 3). Due to 
its presence in TAS and eastern VIC, Rawlinson [29] 
suggested that the delicate skink was a glacial relic 
that had occurred in southern Australia for a pro- 
longed period of time. However, it has been implied 
that the delicate skink might not be native to TAS, 
instead reaching the state via human-assisted colonisa- 
tion. This is because the delicate skink was not 
detected in TAS until 1963, although subsequent 
examination of museum collections revealed that pre- 
viously mis-identified specimens had been collected 
during the 1920s and 1930s [30]. Unlike many reptile 
species whose distribution spans Bass Strait, the deli- 
cate skink does not occur on Wilsons Promontory (the 
most southerly projection of the Australian mainland), 
but it does occur on a Bass Strait Island (Flinders 
Island) that formed part of the Bassian Isthmus during 
the last glacial maxima ([28]; Figure 1, 2). We conduct 
a range of phylogeographic analyses to examine Raw- 
linson's [29] hypothesis, determine the status of the 
Tasmanian population, and investigate the impact of 
historical processes on the evolutionary history of the 
delicate skink. 

Methods 

Sampling 

We obtained tissue samples from 238 Lampropholis 
delicata, representing 120 different populations, from 
across the entire Australian range of the species (Figure 
2, 3; Additional files 1, 2). Samples were obtained from 
the frozen-tissue collections of several Australian 



Museums (Australian Museum, South Australian 
Museum, CSIRO Australian National Wildlife Collec- 
tion), along with our own field collections (Additional 
files 1, 2). We included the closely related L. guichenoti 
(Australian Museum NR2639) and an Australian Eugon- 
.gy/ws-lineage skink Niveoscincus pretiosus (Australian 
Museum NR391) as outgroups in our study. 

DNA extraction, amplification and sequencing 

Total genomic DNA was extracted from liver, muscle, 
toe or tail-tip samples using a Qiagen DNeasy Blood 
and Tissue Extraction Kit (Qiagen, Hilden, Germany). 
For each sample we sequenced portions of four mito- 
chondrial genes: ND2 (-600 bp), ND4 (-700 bp), 
12SrRNA (-700 bp), and 16SrRNA (-500 bp). These 
regions were targeted because work across several taxo- 
nomic levels in squamate reptiles has indicated useful 
levels of variability [e.g. [26,31-34]]. The primers used to 
amplify and sequence these regions are provided in 
Additional file 3. PCR was conducted as outlined in 
Greaves et al. [35], except on a Corbett Research GC1- 
960 thermal cycler. PCR products were purified using 
ExoSAP-IT (USB Corporation, Cleveland, Ohio USA). 
The purified product was sequenced directly using a 
BigDye Terminator v3.1 Cycle Sequencing Kit (Applied 
Biosystems) and then analysed on an ABI 3730XL capil- 
lary sequencer. 

Sequence data were edited using CONTIGEXPRESS in 
VECTOR NTI ADVANCE v9.1.0 (Invitrogen), and 
aligned using the default parameters of CLUSTAL X 
vl.83[36]. We translated all coding region sequences to 
confirm that none contained premature stop codons. 
Sequence data were submitted to GenBank [GenBank: 
JF438009-JF438959, EF567304, EU567726, EU567768, 
EU567769, EU567927, EU567928, EU568019, 
EU568020] (Additional file 1). 

Phylogenetic analyses 

Maximum Likelihood (ML) and Bayesian tree building 
methods were used. We used MODELTEST 3.7 [37] to 
identify the most appropriate model of sequence evolu- 
tion based on the AIC criterion. MODELTEST, con- 
ducted in PAUP* 4.0bl0 [38], was also used to estimate 
base frequencies, substitution rates, the proportion of 
invariable sites (I) and the among-site substitution rate 
variation (G). These values were then used as settings in 
PhyML 3.0 [39] to generate a ML tree with 500 
bootstraps. 

MRBAYES 3.1.2 [40] was then used to complete Baye- 
sian analyses. Preliminary analysis of each mtDNA 
region revealed congruent tree topologies. In order to 
evaluate partitioning strategies, we used MODELTEST 
to determine the most appropriate model for each parti- 
tion. We then conducted a Bayesian analysis for each 
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Figure 2 Map of eastern Australia indicating the collection localities of the Lampropholis delicata samples. The distribution of the nine 
major clades are indicated: clade 1 (grey solid circles), clade 2 (hollow triangles), clade 3 (black solid squares), clade 4 (black solid circles), clade 5 
(grey solid triangles), clade 6 (hollow squares), clade 7 (black solid triangles), clade 8 (hollow circles), clade 9 (grey solid squares). The 
approximate distribution (solid line) of L. delicata is indicated (adapted from [28]), with the population numbers from Additional files 1 and 2 are 
presented next to the sampling localities. 
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Figure 3 The distribution of Lampropholis delicata clades and subclades. The distribution of each subclade identified in Figure 4 and 5 is 
presented (solid lines), along with the location of recognised biogeographic barriers in eastern Australia (dashed lines; Figure 1, Additional file 2). 



partitioning strategy, applying the appropriate model of 
evolution to each partition, and allowing among-parti- 
tion rate variation. We ran each Bayesian analysis for 
five million generations, sampling every 100 generations 
(i.e. 50,000 sampled trees). We ran each analysis twice, 
using four heated chains per run. We discarded the first 



25% of samples as burn-in and the last 37,500 trees 
were used to estimate the Bayesian posterior probabil- 
ities. In order to calculate the AIC and BIC scores for 
different partitions strategies, we calculated the number 
of parameters for each. Following McGuire et al. [41], 
for each parameter we added the number of substitution 
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rates for the model suggested by MODELTEST for that 
partition (6 for GTR, or 2 for HKY), the number of free 
equilibrium base frequencies (3 for GTR and HKY), plus 
one parameter per partition where appropriate for each 
of I and/or G. For multi-partition strategies, we also 
added one parameter per partition, corresponding to the 
among-partition rate multiplier. To calculate the AIC 
and BIC scores, we used the equations: AIC = -2L 2 + Iki 
and BIC = -2L, + (£ 2 ).(ln n) (where i ; is the harmonic 
mean log likelihood for partitioning strategy i, k t is the 
total number of parameters for partitioning strategy i, 
and n is the total number of nucleotides). The program 
TRACER 1.5 [42] was used to check for chain conver- 
gence and mixing. Specifically, raw traces of sampled 
values versus MCMC step numbers were examined to 
confirm that there was no trend away from the mean 
and that there were no large fluctuations in the likeli- 
hood values. 

Bootstrap values (500 ML bootstraps) and Bayesian 
posterior probabilities were used to assess branch sup- 
port. We considered branches supported by bootstrap 
values of 70% or greater [43], and/or posterior probabil- 
ity values greater than or equal to 95% [44] to be sup- 
ported by our data. 

Molecular diversity and population divergence 

Estimates of genetic diversity within L. delicata clades 
(number of haplotypes, h; haplotypic diversity, Hd; num- 
ber of polymorphic sites, S; nucleotide diversity, tt) were 
calculated in DNASP v4.50 [45]. Tamura-Nei (TrN)-cor- 
rected genetic distances within and among clades were 
calculated in MEGA 4 [46]. Genetic differentiation 
among clades within L. delicata was estimated in ARLE- 
QUIN v3.5 [47]. Pairwise O st values (an analogue of 
Wright's fixation index F ST ) were calculated to estimate 
among clade differentiation. We conducted hierarchical 
Analysis of Molecular Variance (AMOVA; [48]) to 
investigate the impact of the a priori (Tables 1 and 2) 
and a posteriori biogeographic barriers on the partition- 
ing of genetic variation within L. delicata. Both tests 
used TrN genetic distances with gamma correction 
(using the value calculated from MODELTEST). Signifi- 
cance levels of all the estimated values were calculated 
by 10,000 permutations, and adjusted according to the 
Bonferroni correction procedure [49] for multiple pair- 
wise comparisons as described by Holm [50]. 

We used Tajima's D [50], Fu's F statistic [51] (calcu- 
lated in ARLEQUIN) and mismatch distributions to test 
for signatures of population expansion within L. delicata 
clades. Significant and negative Tajima's D and Fu's F 
statistic values are indicative of possible population 
expansion. Mismatch frequency histograms were plotted 
in DNASP to determine whether the clades exhibited 
evidence of spatial range expansion or a stationary 



population history [52]. A smooth bell shape signifies 
either population expansion or spatial range expansion, 
whereas a multimodal distribution represents a long his- 
tory in situ [53-56]. To distinguish between these two 
types of distribution, a raggedness index (RI, sum of the 
squared difference between neighbouring peaks) and the 
sum of squared deviations (SSD) between the observed 
and expected mismatch were calculated using the meth- 
ods of Schneider & Excoffier [57] in ARLEQUIN. The 
spatial expansion hypothesis (both RI and SSD) was 
tested using a parametric bootstrap approach (200 
replicates). 

As there are no suitable fossil calibration points 
available for Lampropholis skinks, we estimated the 
divergence time of L. delicata clades using an evolu- 
tionary rate of 1.3-1.63% sequence divergence per mil- 
lion years, based on mitochondrial DNA calibrations 
from other squamate reptile groups (1.3%, [58]; 1.42- 
1.63%, [59]; 1.55%, [60]; 1.62%, [61]; 1.63%, [62]). A 
strict molecular clock (0.0065-0.00815 per lineage sub- 
stitution rate), implemented in BEAST vl.6.1 [63], was 
used to estimate the divergence times within L. deli- 
cata. The Australian Eugongylus lineage is estimated to 
have originated -20 mya [64], and this information 
was used as the maximum age of the tree root. A GTR 
+I+G model of evolution was employed with a coales- 
cent (Bayesian skyline) tree prior. The analysis was run 
twice, with 20 million generations per run (total 40 
million generations). The output was viewed in TRA- 
CER to check that stationarity had been reached, and 
ensure that the effective sample size (ESS) exceeded 
200 [63]. The two separate runs were then combined 
using LOGCOMBINER vl.6.1, with a maximum clade 
credibility tree generated in TREEANNOTATOR 
vl.6.1 and visualised in FIGTREE vl.3.1. A Bayesian 
skyline plot [65] was also generated in TRACER to 
examine the magnitude and timing of population size 
changes in L. delicata. 

Results 

Molecular diversity and phylogeographic structure 

The edited alignment comprised 2426 characters (550 
bp ND2, 671 bp ND4, 708 bp 12SrRNA, 497 bp 
16SrRNA; Additional file 4), of which 813 (33.5%) were 
variable and 587 (24.2%) were parsimony-informative. 
For the ingroup only, the alignment contained 638 
(26.3%) variable characters, of which 543 (22.4%) were 
parsimony-informative. Base frequencies were unequal 
(A = 0.3694, T = 0.2172, C = 0.2889, G = 0.1245), but a 
% 2 test confirmed the homogeneity of base frequencies 
among sequences (df = 498, P = 1.0). The phylogenetic 
analyses were conducted on a dataset comprising the 
165 unique haplotypes that were present within L. deli- 
cata (Table 3). 
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Table 3 Estimates of genetic diversity for the clades present within Lampropholis delicata 



Clade 


n 


h 


Hd 


M(S) 


7T 


Tajima's D 


Fs 


Rl 


SSD 


Clade 1 


14 


13 


0.989 


191(187) 


0.026 


0.378 


0.585 


0.020 


0.033* 


la 


2 


2 


1.0 


38(38) 


0.016 


NA 


NA 


NA 


NA 


lb 


3 


3 


1.0 


7(7) 


0.002 


NA 


NA 


NA 


NA 


1c 


9 


8 


0.972 


19(19) 


0.003 


-0.305 


-2.005 


0.063 


0.054* 


Clade 2 


8 


7 


0.893 


12(12) 


0.001 


-1.576* 


-2.870* 


0.050 


0.010 


Clade 3 


19 


15 


0.965 


170(164) 


0.022 


0.617 


2.306 


0.023 


0.024 


3a 


4 


4 


1.0 


■1('1! 


0.001 


-0.780 


-1.872* 


0.222 


0.042 


3b 


3 


I 


0.0 


0(0) 


0.000 


NA 


NA 


NA 


NA 
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Overall 238 165 0.990 715(625) 0.044 



n = sample size, h = number of haplotypes, Hd = haplotypic diversity, M = total number of mutations, S = number of segregating (polymorphic) sites, n = 
nucleotide diversity, Fs = Fu's F statistic, Rl = raggedness index, SSD = sum of squared deviations. Asterisks indicate significant Tajima's D, Fu's F statistic, Rl and 
SSD values. 



The AIC from MODELTEST supported the GTR + I + 
G substitution model as the most appropriate for our 
unpartitioned dataset. Parameters estimated under this 
model were: relative substitution rates (A<->C = 2.1053, 
A<->G = 44.1898, A^T = 2.5350, C^G = 1.1020, C^T 
= 29.2135, relative to G<->T = 1.00), proportion of 
invariable sites (0.5241), and gamma distribution shape 
parameter (0.7015). We evaluated three partitioning 
strategies for our dataset (Table 4). The unpartitioned 
and by gene partitioning strategy were analysed using 
the GTR + I + G model for all nucletotides. When the 
data was partitioned by codon, we used a mixture of 
GTR and HKY models for each partition (Table 4). Both 
the AIC and BIC scores ranked the most highly parame- 
terised strategy (by gene and codon) as the most appro- 
priate. However, the topologies of the ML, unpartitioned 
Bayesian and partitioned Bayesian trees were congruent, 
therefore we present the optimal ML tree (-In L = 
14501.43296) with ML bootstrap (BS) values and unpar- 
titioned Bayesian posterior probabilities (PP) indicating 



branch support (Figure 4, 5). Nine well-supported, non- 
overlapping clades (labelled from the top to bottom of 
the tree, and roughly related to their north to south 
geographic distribution) are present within L. delicata 
(Figure 2, 3, 4, 5), with high levels of haplotypic and 
nucleotide diversity within each clade (Table 3). The 
PP's of the main clades and subclades in the partitioned 
Bayesian analysis were identical to those from the 
unpartitioned Bayesian analysis presented in Figure 4 
and 5, except that the support value for Clade 9 was 
lower (0.85 rather than 0.95). 

Clade 1 encompasses populations from coastal north- 
ern and eastern QLD, and is comprised of three sub- 
clades (Figure 2, 3, 4, 5). Subclade la includes 
populations north of Townsville, subclade lb contains 
populations in the Mackay region, and subclade lc 
stretches from the Rockhampton region to Bania State 
Forest (inland from Bundaberg) in south-east QLD (Fig- 
ure 2, 3, 4, 5). Clade 2 is restricted to Kroombit Tops 
(Figure 2, 3, 4, 5). A complex mosaic of geographically 
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Table 4 Test of alternative partitioning strategies 

Partition Strategy No. of Partitions No. of parameters (k,) Harmonic Mean Log-likelihood (/.,) AIC Score BIC Score 

Unpardoned 1 11 -15185.35 30392.7 30544.7 

By gene 4 48 -14753.32 29602.6 30266.0 

By gene and codon 8 82 -14448.71 29061.4* 30194.6* 

Modeltest selected the GTR+I+G model for all partitions except for codon position 1 in ND2 and ND4 (HKY+G), ND2 codon position 2 (GTR+I), ND2 codon 
position 3 and ND2 codon positions 2 & 3 (GTR+G). The asterisk indicates the optimal partitioning strategy. 



non-overlapping clades occurs throughout south-eastern 
QLD and northern NSW (Figure 2, 3, 4, 5). Clade 3 
includes populations from the Sunshine Coast and Bris- 
bane region of south-eastern QLD (Figure 2, 3, 4, 5). 
There is strong support for a close affinity between 
Clade 2 and Clade 3 (Figure 4, 5). Subclade 3a occurs 
within Cooloola National Park, while subclade 3b is 
restricted to the Bunya Mountains (Figure 2, 3, 4, 5). 
Subclade 3c extends from the Maryborough region to 
the northern suburbs of Brisbane, with subclade 3d con- 
taining populations from the southern suburbs of Bris- 
bane to Barney View on the western side of Lamington 
National Park (Figure 2, 3, 4, 5). 

Clade 4 (96 BS) occurs along the Main Range and 
throughout the QLD/NSW border region, and in inland 
northern NSW (Figure 2, 3, 4, 5). Subclade 4a (93 BS) 
extends from Deongwar State Forest and other areas in 
south-eastern QLD (Lamington NP, Main Range) 
through the elevated regions of inland northern NSW to 
Riamukka State Forest, inland from Port Macquarie 
(Figure 2, 3, 4, 5). Subclade 4b is restricted to the more 
coastal regions of northern NSW (Mt Warning NP, Bor- 
der Ranges NP, Nightcap NP, Whian Whian State For- 
est, Alstonville), while subclade 4c occurs further inland 
at Bolivia Hill, Torrington State Forest and near Armi- 
dale (Figure 2, 3, 4, 5). Clade 5 extends along the north- 
ern NSW coastal region from Yamba to Cairncross 
State Forest near Port Macquarie, and is comprised of 
two subclades: subclade 5a (Yamba to Coffs Harbour) 
and subclade 5b (Styx River State Forest to Cairncross 
State Forest) (Figure 2, 3, 4, 5). Clade 6 represents a dis- 
junct inland population on the Blackdown Tableland in 
southern QLD (Figure 2, 3, 4, 5). 

Clade 7 is geographically widespread, occurring from 
the Australian Capital Territory (ACT) and inland 
southern NSW (subclade 7a) across to western VIC (Lit- 
tle Desert NP) and south-eastern SA (subclade 7b) (Fig- 
ure 2, 3, 4, 5). Clade 8 is a disjunct population that 
occurs in Coolah Tops National Park in inland northern 
NSW (Figure 2, 3, 4, 5). Clade 9 is distributed from the 
central coast of NSW and throughout eastern VIC and 
TAS (Figure 2, 3, 4, 5). Subclade 9a encompasses the 
central coast of NSW and the Sydney region, subclade 
9b occurs at Brayton, while subclade 9c comprises 
populations from coastal southern NSW (Figure 2, 3, 4, 



5). Subclade 9d represents a shallow clade that is dis- 
tributed throughout eastern VIC and TAS (Figure 2, 3, 
4, 5). 

Genetic differentiation among clades and divergence time 
estimates 

Considerable genetic differentiation was evident 
amongst the nine L. delicata clades, with extremely high 
and statistically significant pairwise O st values among 
clades (Table 5). The only comparisons that were not 
significant were those involving clades with low sample 
sizes (e.g. Clade 8). Substantial genetic distances are evi- 
dent among the clades (4.3-7.4%; Table 5), indicating 
that the divergences within L. delicata occurred during 
the late Miocene-Pliocene (Figure 6, Additional file 5). 
The intra-clade genetic divergences in L. delicata were 
0.0-2.6% (Table 3). The vast majority (97.7%) of genetic 
variation in L. delicata was partitioned among popula- 
tions (Table 6). The nine a priori biogeographic barriers 
(Tables 1 and 2) accounted for 64.1% of the genetic var- 
iation in L. delicata (Table 6). This value increased to 
66.5% when the two barriers identified a posteriori 
(Blackdown Tableland, Coolah Tops) were included in 
the analysis (Table 6). 

The Bayesian skyline plot indicated recent (i.e. last 0.2 
myr) contraction then expansion of L. delicata popula- 
tions (Figure 7), although there was no consistent sup- 
port for the model of spatial expansion in L. delicata 
clades or subclades. Three main clades (2, 4 and 7) and 
four subclades (3a, 4b, 7b and 9d) deviated significantly 
from the expectations of neutrality (Tajima's D, Fu's F 
statistic) (Table 3), suggesting recent population expan- 
sion. However, the RI and SSD values indicated that a 
model of population expansion could only be conclu- 
sively rejected for two main clades (1 and 4) and two 
subclades (lc and 4c) (Table 3). 

Discussion 

The nine main clades of Lampropholis delicata appear 
to have diverged during the late Miocene-Pliocene. 
Although the current study relied solely on mitochon- 
drial DNA sequence data, the same tree topology is evi- 
dent in a molecular phylogeny for the Lampropholis 
genus (including representatives from each main L. deli- 
cata clade) based on mitochondrial DNA and five 
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Figure 4 Phylogram for the delicate skink [Lampropholis delicata). The phytogeny is based on 2426 bp of mitochondrial DNA (550 bp ND2, 
671 bp ND4, 708 bp USrRNA, 497 bp WSrRNA). The population numbers (Figure 2, Additional file 2) are provided in parentheses. The overall 
tree topology is indicated in the inset. Nine major genetic clades are identified within the delicate skink. Two measures of branch support are 
indicated with ML bootstraps (500 replicates) on the left and Bayesian posterior probabilities on the right (only values over 70 and 0.9, 
respectively, are shown). 
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Figure 5 Continuation of Figure 4. 



nuclear genes (C. Hoskin, C. Moritz & D. Chappie, 
unpublished data). The divergence of L, delicata corre- 
sponds to a time when rainforest habitat in eastern Aus- 
tralia was in decline as a result of a drying climate, 
resulting in restriction of rainforest to a series of dis- 
junct remnants that have been described as an 'archipe- 
lago of refugia' [1,3,17]. The delicate skink occurs in 
rainforest or rainforest fringes and therefore likely 



experienced similar reduction and fragmentation, result- 
ing in genetic divergence among geographically isolated 
populations. Despite evidence for the expansion and 
contraction of some clades throughout the Pleistocene, 
each is geographically structured and non-overlapping 
(Figure 2, 3). This pattern that has been observed in a 
range of other taxa (Table 2; [1]), including L. guichenoti 
[26]. Phylogeographic breaks in the delicate skink 
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Table 5 Mean Tamura-Nei corrected mtDNA genetic 
distances (below diagonal) and pairwise $ ST (above 
diagonal) among the major clades (1-9) identified in 
Figure 4 and 5 
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0.902* 
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0.060 


0.055 
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9 


0.070 


0.053 


0.050 


0.050 


0.050 


0.050 


0.054 


0.049 





Asterisks denote statistical significance following Bonferroni correction. 
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Figure 6 BEAST maximum clade credibility tree for 
Lampropholis delicata The divergence times correspond to the 
mean posterior estimate of their age in millions of years. The bars 
indicate the 95% HPD interval for the divergence time estimates. 
Nodes with a posterior probability of > 0.95 are indicated by an 
asterisk. The clades and subclades identified in Figure 4 and 5 are 
presented. 



generally correspond to dry habitat and topographic bar- 
riers (i.e. Burdekin Gap, St Lawrence Gap, Kroombit 
Tops, McPherson Range; Hunter Valley, southern NSW; 
Figure 1, 2, 3). However, contrary to the hypothesis of 
Rawlinson [29], the delicate skink appears to be a rela- 
tively recent arrival in south-eastern Australia and exhi- 
bits no evidence of restricted geneflow across the 
barriers in this region (e.g. East Gippsland, Bass Strait). 

Phylogeographic structure in the delicate skink 
corresponds to dry habitat and elevational barriers 

Despite its widespread distribution along the east coast 
of Australia, there is substantial phylogeographic struc- 
ture across the native range of the delicate skink (Figure 
4, 5). In many instances these breaks are concordant 
with dry habitat corridors, indicating that regions of 
drier vegetation represent effective barriers to dispersal 
for the mesic-adapted delicate skink. For instance, the 
delicate skink exhibits a moderate genetic break (4.5%, 
subclade la vs lb; Figure 6) across the Burdekin Gap in 
North QLD. Equivalent Pliocene divergences between 
populations either side of the Burdekin Gap have been 
reported in open forest frogs [66,67], rainforest lizards 
[22], woodland lizards [68], rainforest birds (Pleistocene; 
[19]) and freshwater fish (Miocene; [69]) (Table 2). In 
contrast, the St Lawrence Gap north of Rockhampton 
on the central QLD coast has only been identified as a 
significant barrier for one lizard species (late Pleisto- 
cene-early Pliocene, [68]; Table 2). Although divergence 
is also evident across the St Lawrence Gap in the deli- 
cate skink (4.8%, mid-late Pliocene, subclade lb vs lc), a 
more substantial break is evident a little to the south, 
between clades 1 and 3 (7.1%, late Miocene) in the 
Gladstone region (Table 5, Figure 2, 3, 4, 5). 

Two high elevation areas (Kroombit Tops, Blackdown 
Tableland) inland from the main range of L. delicata in 
southern QLD were found to harbour genetically diver- 
gent lineages. Both areas are remnant patches of moist 
forest that are surrounded by drier lowland eucalypt 
woodland, and are disjunct from the main distribution 
of the delicate skink along the east coast ([70]; Figure 2, 
3). Kroombit Tops (-730 m) was identified a priori as a 
potential habitat isolate for the delicate skink, as the 
region provides a cooler and wetter refuge for mesic- 
adapted species ([70,71]; Tables 1 and 2). The Kroombit 
Tops population of the delicate skink diverged from the 
surrounding coastal populations during the mid Pliocene 
(4.9%, Clade 2 vs 3; Table 5, Figure 6), a pattern that has 
also been observed in a rainforest bird [71], two open 
forest frogs [66,67], and several open forest reptiles [72] 
(Table 2). The Blackdown Tablelands are a moderate 
elevation plateau (-600 m) that provides an isolated 
refugium for numerous mesic-adapted species. Our 
results for the delicate skink (4.3-6.7%, Pliocene; Figure 
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Table 6 Hierarchical Analysis of Molecular Variance (AMOVA) for Lampropholis delicata populations and biogeographic 
regions 



Comparison 




Observed Partition 




Total Variance 
Components 


Total Sum of 
Squares 


P 




Among 
Populations 


Within Populations 










All populations 


97.7% (118) 

Among 
Regions 


2.3% (1 1 9) 

Among Populations 
within Regions 


Within 
Populations 


70.09 


1 6488.29 


< 0.001 (237) 


a priori biogeographic 
barriers 


64.1% (9) 


33.8% (109) 


2.1% (119) 


78.72 


1 6488.29 


All comparisons < 
0.001 (237) 


a priori + a posteriori 
biogeographic barriers 


66.5% (11) 


31.4% (107) 


2.1% (119) 


78.57 


1 6488.29 


All comparisons < 
0.001 (237) 



Separate analyses are conducted for the nine biogeographic barriers that were identified a priori (Tables 1 and 2), and including the two additional barriers 
identified a posteriori (i.e. Blackdown Tableland, Coolah Tops). The degrees of freedom are indicated in parentheses. Statistical significance (P) was tested with 
10,000 permutations. 



6) provide evidence that the fauna of this region may 
also be genetically divergent from the coastal popula- 
tions, a pattern also seen in other open forest reptiles 
[72]. 

The Coolah Tops are a high elevation (1200 m) pla- 
teau located in inland northern NSW, just to the north 
of the Hunter Valley (Figure 2, 3). The refugial popula- 
tion of the delicate skink that occurs on the Coolah 
Tops was found to have diverged from the nearby popu- 
lations in northern NSW during the early-mid Pliocene 
(4.6-4.8%; Table 5, Figure 6). The dry habitat in the 
Hunter River Valley has been demonstrated to represent 
a major barrier to dispersal in both woodland and wet 
forest species ([73-76]; Table 2). While the divergence 
across the Hunter Valley was estimated to have 
occurred in the Miocene in the congeneric L. guichenoti, 
which is frequently sympatric with L. delicata [26], an 
early-mid Pliocene break was observed across this bar- 
rier in the delicate skink (subclade 9a vs Clades 4-5; Fig- 
ure 2, 3, 4, 5). The divergence estimate for the delicate 




Time (million years ago) 

Figure 7 Bayesian skyline plot for Lampropholis delicata. The 

estimated effective population size over the past 6 myr is 
presented. The solid line indicates the median population estimate, 
with the two light blue lines indicating the 95% confidence interval. 



skink is consistent with those reported for most other 
species across the Hunter Valley (Table 2). 

Our analyses revealed a complex mosaic of geographi- 
cally structured, non-overlapping clades and subclades 
(Clades 3-5) in the delicate skink in south-eastern QLD 
and northern NSW (Figure 2, 3, 4, 5). The McPherson 
Range that occurs along the border region (Figure 1) is 
concordant with the phylogeographic break between 
Clades 3 and 4 (Figure 2, 3). The distribution of Clade 4 
extends northwards into south-east QLD to the Main 
Range, which runs perpendicular to the western edge of 
the McPherson Range (Figure 1). Similar biogeographic 
patterns involving the Main and McPherson Ranges 
occur in wet forest {Litoria pearsoniana; [77]) and open 
forest frogs {Litoria fallax; [66]) (Table 2). The early- 
mid Pliocene split found across the McPherson/Main 
Ranges in the delicate skink (4.8%; Figure 6) is concor- 
dant with that observed in L. guichenoti, which also 
inhabits open woodlands and dry sclerophyll forest [28], 
but intermediate between that reported for frogs (Mio- 
cene; [66,77]) and a wet forest snake (Pleistocene; [78]) 
(Table 2). 

Some relatively minor phylogeographic structure is 
evident among the populations in the Maryborough, 
Sunshine Coast and Brisbane regions of south-east QLD 
(subclades 3a-d; Figure 2, 3, 4, 5). A more substantial 
break occurs between inland (Clade 4) and coastal 
(Clade 5) delicate skink populations north of the Hunter 
Valley in NSW (4.5%, early-mid Pliocene; Table 5, Fig- 
ure 6). A similar coastal vs inland divergence in north- 
ern NSW is shared with White's skink (Liopholis whitii, 
[74]), and reflects an equivalent pattern that is regularly 
observed in southern NSW (Table 2). Indeed, an analo- 
gous pattern is evident within Clade 4 in northern 
NSW, with subclade 4b occurring near the coastal mar- 
gin, subclade 4a present in intermediate areas (with sec- 
ondary contact between 4a and 4b occurring in the 
Border Ranges NP), and subclade 4c occurring further 
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inland (Figure 2, 3). These biogeographic patterns, com- 
bined with the break observed in southern NSW (Figure 
2, 3), indicate that high elevation areas may represent 
barriers to dispersal in the delicate skink. 

The delicate skink is a relatively recent arrival in southern 
Australia 

The five phylogeographic studies that have had sufficient 
sampling to examine the impact of the elevational and 
habitat barriers in southern NSW have reported a 
genetic break between the inland (including the ACT) 
and coastal regions ([26,71,73,74,79]; Table 2). The 
impact of this barrier is pronounced in the delicate 
skink, as populations from the ACT and inland southern 
NSW are more closely related to the SA populations 
(3.0%, early Pleistocene-late Pliocene) than the adjacent 
populations along the NSW coast (5.3%, mid-Pliocene; 
Figure 2, 3, 4, 5). This indicates that the delicate skink 
most likely reached SA from the southern NSW region 
via an inland route, rather than along a coastal dispersal 
pathway through VIC. The delicate skink may have dis- 
persed through the mesic vegetation that is located 
along the Murray River, which forms the border 
between NSW and VIC for the majority of its length 
(Figure 2, 3). Indeed, the eastern water skink (Eulam- 
prus quoyii) has a continuous distribution through the 
Murray-Darling River system that connects populations 
along the east coast (North QLD to southern NSW) to 
an isolated population in south-eastern SA [28]. This 
biogeographic pattern was not previously suspected for 
the delicate skink and explains the large distributional 
gap across western VIC between the eastern suburbs of 
Melbourne and Little Desert NP in north-western VIC 
(Figure 2, 3). Given this pattern, it was not possible to 
examine the impact of the Murray Basin on the delicate 
skink (Table 2). 

Several frog and lizard species exhibit deep phylogeo- 
graphic breaks in the East Gippsland region 
[26,67,74,79], a pattern that is believed to be the result 
of repeated marine inundation of the area since the 
Miocene (Tables 1 and 2). However, the East Gippsland 
region does not appear to constitute a barrier to disper- 
sal in the delicate skink, with an extremely shallow clade 
(intraclade divergence 0.2%) distributed across eastern 
VIC and across Bass Strait into TAS (Figure 2, 3). The 
coastal area in East Gippsland has been relatively stable 
since the late Pleistocene [15,80], enabling the delicate 
to disperse across eastern VIC from southern NSW. 

The delicate skink has colonised TAS during the late 
Pleistocene, with the presence of shared haplotypes 
between populations in eastern VIC and TAS indicating 
a connection between these two regions until relatively 
recently (~12-15 kya; Table 5, Figure 6, Additional file 
2). The timing coincides with the inundation of the 



Bassian Isthmus, which connected eastern VIC (Wilsons 
Promotory) to north-eastern TAS between 43-14 kya 
([16]; Table 1, Figure 1). Although the delicate skink 
does not currently occur on Wilsons Promontory, it is 
present on Flinders Island (which formed part of the 
Bassian Isthmus) and is common throughout north-east- 
ern TAS ([[28,81], DGC, personal observation]). While a 
second land bridge was located from western VIC, 
through the King Island region to western TAS from 
~43-17.5 kya [16], the absence of the delicate skink 
from western VIC would have precluded dispersal of the 
species via this western route. Fossil evidence for a 
Nothofagus tree species on King Island 38 kya suggests 
that moist forest vegetation occurred along the Bass 
Strait land bridges [82], enabling dispersal of the delicate 
skink into TAS. Although some other species appear not 
to have used these recent land bridges (frog: Crinia sign- 
igera, [79]; reptiles: Acritoscincus duperreyi, [76]; Lerista 
bougainvilli, [83]; mammals: Dasyurus maculatus, [84]), 
others appear to have dispersed across these Bass Strait 
land bridges (frogs: Limnodynastes peronii and tasma- 
niensis, [67]; reptiles: Liopholis whitii, [73]; Notechis scu- 
tatus, [85]) (Table 2). The repeated presence of the land 
bridges has also restricted east-west gene flow across 
Bass Strait in several marine invertebrate species 
([86-89]; Table 2). 

There is no evidence to support the hypothesis of 
Rawlinson [29] that the delicate skink is a 'glacial relic' 
with a relatively long presence in southern Australia. In 
contrast, our analyses indicate that the delicate skink 
only colonised VIC and TAS during the late Pleistocene 
from coastal southern NSW. Although the delicate 
skink (also known as the rainbow or plague skink in its 
introduced range) is a successful invasive species in the 
Hawaiian Islands, New Zealand, and Lord Howe Island 
[90], there is no strong evidence to suggest that it repre- 
sents an introduced species in TAS. However, given the 
relatively shallow genetic divergences within subclade 
9d, we are unable to completely exclude the possibility 
that the delicate skink reached TAS via human-asso- 
ciated colonisation. 

Conclusions 

We performed a detailed phylogeographic study of a 
species found in mesic forests down almost the entire 
length of eastern Australia. Lampropholis delicata is a 
single widespread, but genetically variable, species con- 
sisting of nine geographically structured, non-overlap- 
ping clades. This structuring is likely the result of 
population subdivision across dry habitat barriers (Bur- 
dekin Gap, St Lawrence Gap, Hunter Valley), topo- 
graphic barriers (McPherson Range, southern NSW) and 
to upland habitat isolates (Kroombit Tops, Blackdown 
Tableland, Coolah Tops). In contrast, in the south-east 
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of its range, the delicate skink exhibits evidence for 
recent dispersal into SA via an inland route, and 
through eastern VIC and across the Bassian Isthmus 
into TAS. Previous studies have demonstrated geo- 
graphic variation in morphology, reproductive ecology 
and life history in the delicate skink [91,92]. Given the 
presence of multiple divergent lineages across the range, 
this regional variation in morphology and life history 
may have a genetic, as well as climatic or environmental, 
basis. 

Additional material 



Additional file 1: Complete collection locality table, with museum 
specimen and tissue voucher number and GenBank accession 
numbers 

Additional file 2: Clades, haplotypes, latitude and longitude for 
Lampropholis delicata populations sampled in the study 

Additional file 3: Oligonucleotide primers used in this study 

Additional file 4: The concatenated alignment (fasta format) for the 
165 Lampropholis delicata haplotypes (ND2: 1-550, ND4: 551-1221, 
12SrRNA: 1222-1929, 16SrRNA: 1930-2426) 

Additional file 5: Divergence time estimates for the main 
Lampropholis delicata clades and subclades 
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